This is file 4 of an exploratory data analysis (EDA) of the NEFSC Bottom Trawl Survey catch counts from 2008 through 2019.

EDA 1 cleans the species data using FishBase and SeaLifeBase.
EDA 2 lists out all the species names, does some more name cleaning, and combines the species information with the trawl count data. EDA 3 queries FishBase for swim bladder information (but I do not trust it) and SeaLifeBase and FishBase for habitat information.

This EDA plots out the spatial and temporal study areas, finds a few more data mistakes, and takes a high-level look at the counts of Fish and Not Fish species.

knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = FALSE, 
                      fig.width = 7, fig.height = 6.5)
library(knitr)
library(tidyverse)
library(lubridate)
library(kableExtra)
library(rfishbase)
library(sf)
library(leaflet)

Loading the data

These are the outputs from the first and third EDA files.

dta <- read_rds("output_input/dta.rds")
dta_ed <- read_rds("output_input/dta_ed.rds")

sp_fish <- read_rds("output_input/sp_fish_2.rds") 
sp_sealife <- read_rds("output_input/sp_sealife_2.rds")

non_sp <- read_rds("output_input/non_sp.rds")
non_sp_habitat <- read_rds("output_input/non_sp_habitat.rds")

Spatial data

Creating a spatial data frame to plot the stations.

meta_data <- c("cruiseid", "cruise", "station", "beginlon", "beginlat", "begindatetime", "stratum")

dta_sp <- dta %>%
  filter(!is.na(begindatetime)) %>% 
  select(all_of(meta_data)) %>% 
  st_as_sf(coords = c("beginlon", "beginlat"), crs = 4326)

strata <- st_read("shpfiles/BTS_strata/no_prj/finstr.shp", crs = 4326)
## Reading layer `finstr' from data source `C:\Users\ealab\Documents\research_projects\TFO\R_bot_trawl_survey\shpfiles\BTS_strata\no_prj\finstr.shp' using driver `ESRI Shapefile'
## Simple feature collection with 310 features and 6 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: -79.18422 ymin: 32.50221 xmax: -65.16869 ymax: 44.83202
## geographic CRS: WGS 84

Basemap creation

base_map <- leaflet(options = leafletOptions(zoomControl = TRUE,
                                             zoomSnap = 0.1)) %>% 
  # add ocean basemap
  addProviderTiles("Esri.OceanBasemap") %>%

  # add another layer with place names
  addProviderTiles("Hydda.RoadsAndLabels", group = 'Place names') %>%
  
  # add graticules from a NOAA webserver - commented out because I think they
  # add clutter
  # addWMSTiles(
  #   "https://gis.ngdc.noaa.gov/arcgis/services/graticule/MapServer/WMSServer/",
  #   layers = c("1-degree grid", "5-degree grid"),
  #   options = WMSTileOptions(format = "image/png8", transparent = TRUE),
  #   attribution = NULL,group = 'Graticules') %>%
  
  # add a map scalebar
  addScaleBar(position = 'bottomright') %>% 

# focus map in a certain area / zoom level
  setView(lng = -72.201, lat = 39.308, zoom = 6)

# or set the boundary
#fitBounds(lng1 = -82.4, lat1 = 16.0, lng2 = -32.0, lat2 = 71.2)

Stations overlaid on strata

The dots have a degree of transparency - darker blue areas indicate more sampling than lighter blue areas.

base_map %>% 
  addPolygons(data = strata,
              color = "#000000",
              weight = 1,
              fillColor = "C0C0C0") %>% 
  addCircles(data = dta_sp,
             weight = 2) %>% 
  addMiniMap(position = "topleft", width = 125, height = 125)


I removed the projection file (world Mercator) from the shapefile before reading it in. This made it a lot easier to give the strata polygon a coordinate reference system (CRS) that was not tainted by a projection.

Quick check that the point data and polygon data match up.

I was playing with the data in QGIS and found some discrepancies between strata ID and strata listed in the BTS data file. I initially thought the strata discrepancies had to do with the World Mercator projection (BTS file is not projected), but after removing the projection and making sure both the point data and polygon data have the same CRS, there are still strata discrepancies. The following is a list of BTS points with a stratum that does not match the stratum of the polygon it is in.

# pulling out the attributes from the polygon to add to the points
strata_id <- strata %>% 
  select(FINSTR_G_, FINSTR_G_I, STRATA)

# Spatial Intersection 
# adding the polygon information to the points
id_check <- dta_sp %>% 
  st_intersection(strata_id) %>% 
  filter(stratum != STRATA)

# printing out the points with strata information that does NOT match 
# the strata id in the polygon file
id_check %>%
  kable(caption = "These data points from the BTS data sets might have mislabeled strata information. Blue = BTS data. Grey = information sampled from strata polygons.") %>% 
  kable_styling(bootstrap_options = c("striped", "hover")) %>% 
  column_spec(c(5, 8), bold = TRUE) %>% 
  column_spec(1:5, color = "blue") %>% 
  column_spec(6:9, color = "grey")
These data points from the BTS data sets might have mislabeled strata information. Blue = BTS data. Grey = information sampled from strata polygons.
cruiseid cruise station begindatetime stratum FINSTR_G_ FINSTR_G_I STRATA geometry
HB201702 201702 344 2017-05-08 01:41:46 1340 15 1351 1351 POINT (-66.81511 44.0805)
HB201002 201001 341 2010-04-22 06:23:20 1400 45 3750 3750 POINT (-70.10708 43.49685)
HB201601 201602 371 2016-06-03 03:48:27 3660 51 1260 1260 POINT (-70.48854 42.61377)
HB201401 201402 331 2014-05-28 19:15:04 3650 62 3640 3640 POINT (-70.61683 42.25367)
HB201601 201602 314 2016-05-27 19:30:55 1220 66 1290 1290 POINT (-66.75437 42.21039)
HB200802 200804 224 2008-04-09 10:04:27 1220 73 1210 1210 POINT (-66.29068 42.1382)
HB201702 201702 316 2017-05-03 12:52:52 3560 76 3610 3610 POINT (-70.084 42.08083)
HB201002 201001 383 2010-04-30 08:39:56 3590 78 3600 3600 POINT (-70.15049 41.95197)
HB200901 200902 288 2009-04-17 21:57:18 1180 87 1170 1170 POINT (-66.43684 41.03883)
HB201301 201302 274 2013-04-17 10:04:02 1180 87 1170 1170 POINT (-65.73796 41.77114)
HB201002 201001 281 2010-04-15 12:42:52 1200 91 1230 1230 POINT (-68.47917 41.28285)
HB201002 201001 195 2010-03-24 17:26:21 3020 109 3450 3450 POINT (-71.59298 41.08289)
HB200901 200902 171 2009-03-24 21:54:52 3450 111 1050 1050 POINT (-71.2306 41.34053)
HB200901 200902 166 2009-03-24 12:10:24 3020 111 1050 1050 POINT (-71.71487 41.03952)
HB201102 201102 710 2011-03-31 04:46:56 1500 116 1090 1090 POINT (-70.55168 40.91987)
HB201102 201102 709 2011-03-31 04:26:01 1500 116 1090 1090 POINT (-70.5396 40.90333)
HB200901 200902 1 2009-02-27 23:08:00 3050 120 3020 3020 POINT (-71.98233 40.9425)
HB200901 200902 4 2009-02-28 06:42:58 3050 120 3020 3020 POINT (-71.99636 40.93274)
HB201102 201102 136 2011-03-22 01:21:44 3050 120 3020 3020 POINT (-72.00169 40.93943)
HB201802 201802 78 2018-03-30 10:02:56 3170 140 3140 3140 POINT (-73.8403 40.17061)
HB201401 201402 166 2014-04-22 01:00:53 1120 144 1110 1110 POINT (-69.26267 39.9895)
HB201102 201102 176 2011-03-27 13:19:14 1120 154 0 0 POINT (-70.36648 39.96355)
HB201901 201902 94 2019-03-19 10:23:36 1750 155 1030 1030 POINT (-72.74958 39.21102)
HB201002 201001 149 2010-03-20 09:51:02 1750 179 1760 1760 POINT (-72.96657 38.8439)
HB201002 201001 150 2010-03-20 11:41:28 1750 179 1760 1760 POINT (-72.96685 38.83938)
HB201102 201102 104 2011-03-13 20:38:22 1750 179 1760 1760 POINT (-72.96497 38.85372)
HB201401 201402 37 2014-04-01 04:43:35 3290 180 1690 1690 POINT (-74.9673 38.00234)
HB201002 201001 102 2010-03-10 11:22:48 1720 186 1710 1710 POINT (-73.72298 38.1579)
HB201102 201102 2 2011-03-03 16:07:10 1720 186 1710 1710 POINT (-73.57182 38.36966)
HB201501 201502 125 2015-03-28 08:29:49 1720 186 1710 1710 POINT (-73.37072 38.44075)
HB201802 201802 55 2018-03-20 06:25:38 1720 186 1710 1710 POINT (-73.86436 38.04671)
HB201901 201902 77 2019-03-17 21:56:03 1720 186 1710 1710 POINT (-73.57659 38.36939)
HB201201 201202 112 2012-03-22 17:47:17 3260 189 3290 3290 POINT (-74.87447 38.49778)
HB201002 201001 146 2010-03-20 01:25:05 1710 190 1720 1720 POINT (-73.47284 38.39103)
HB201002 201001 84 2010-03-08 20:10:49 1670 202 1660 1660 POINT (-74.65666 37.02749)
HB201301 201302 22 2013-03-17 14:05:57 1680 203 1670 1670 POINT (-74.70588 36.63426)
HB201601 201602 60 2016-04-12 15:54:00 1680 203 1670 1670 POINT (-74.64055 36.86818)
HB201802 201802 44 2018-03-18 22:49:49 1680 203 1670 1670 POINT (-74.66891 36.71928)
HB201802 201802 47 2018-03-19 09:22:55 1680 203 1670 1670 POINT (-74.71052 37.10733)
HB201901 201902 71 2019-03-17 07:09:40 1720 205 1680 1680 POINT (-74.32036 37.53059)
HB201802 201802 37 2018-03-18 06:36:59 1640 230 1620 1620 POINT (-74.7835 36.09291)
HB200901 200902 36 2009-03-07 15:49:11 1620 231 1630 1630 POINT (-74.85443 35.66398)
HB200901 200902 33 2009-03-07 09:21:14 1640 231 1630 1630 POINT (-74.77941 36.10063)
HB200901 200902 32 2009-03-07 06:44:14 1640 231 1630 1630 POINT (-74.76769 36.14999)
HB201002 201001 75 2010-03-07 22:04:01 1620 231 1630 1630 POINT (-74.80907 36.42635)
HB201002 201001 70 2010-03-07 09:57:05 1640 231 1630 1630 POINT (-74.78203 36.05996)
HB201201 201202 51 2012-03-05 19:39:10 1640 231 1630 1630 POINT (-74.78047 36.06017)
HB201201 201202 52 2012-03-05 21:05:30 1640 231 1630 1630 POINT (-74.7826 36.05915)
HB200901 200902 37 2009-03-07 17:48:06 1620 232 1640 1640 POINT (-74.84881 35.51938)
HB201102 201102 23 2011-03-05 23:22:28 1630 232 1640 1640 POINT (-74.7856 36.39567)
HB201301 201302 53 2013-03-20 10:01:27 1630 232 1640 1640 POINT (-74.81488 35.65173)
HB201501 201502 72 2015-03-20 20:09:37 1630 232 1640 1640 POINT (-74.78115 36.42189)
HB201702 201702 31 2017-03-10 02:22:23 1630 232 1640 1640 POINT (-74.7721 36.45907)
HB201702 201702 34 2017-03-10 07:12:04 1630 232 1640 1640 POINT (-74.78428 36.39262)
HB201802 201802 39 2018-03-18 12:20:40 1630 232 1640 1640 POINT (-74.79264 36.26113)
HB200901 200902 60 2009-03-09 13:36:23 3440 237 3430 3430 POINT (-75.41798 35.77488)
HB201601 201602 71 2016-04-13 17:39:27 3410 238 3440 3440 POINT (-75.32338 35.82413)
HB201601 201602 83 2016-04-14 21:05:29 8520 243 8999 8510 POINT (-75.54844 34.76052)
HB201102 201102 48 2011-03-09 02:40:51 8500 244 8999 8520 POINT (-74.87597 35.4742)
HB201301 201302 48 2013-03-20 00:39:16 8520 245 8999 8530 POINT (-74.89261 35.40634)
HB201002 201001 54 2010-03-05 17:42:17 7520 265 755 7550 POINT (-75.91214 35.01254)


When viewing map below in your browser, you can zoom in and out. I also added the point attributes to a popup. Click on a station and the data will pop up.

base_map %>% 
   addPolygons(data = strata,
              color = "#000000",
              weight = 1,
              fillColor = "C0C0C0") %>% 
  addCircles(data = id_check,
             weight = 4,
             opacity = 1,
             popup = ~paste0("Cruise ID: ", cruiseid,"<br/>", 
                             "Station: ", station,"<br/>", 
                             "Datetime: ", begindatetime, "<br/>", 
                             "BTS stratum: ", stratum, "<br/>",
                             "Polygon stratum: ", STRATA, "<br/>",
                             "long/lat: ", geometry)) %>%
  addMiniMap(position = "topleft", width = 125, height = 125)


Questions:

  • How were the strata determined?
    • Are they an ecological unit based on depth?
  • Should the strata for the above BTS points be changed to the corresponding the polygon strata?

Temporal data

Each dot is the datetime and yearday of a sampling period (station).

dta_1 <- dta %>%
  filter(!is.na(begindatetime)) %>% 
  select(all_of(meta_data)) %>% 
  mutate(yearday = yday(begindatetime),
         month = month(begindatetime, label = TRUE, abbr = FALSE),
         year = year(begindatetime))

start <- dta_1 %>% 
  filter(!is.na(begindatetime)) %>% 
  summarize(min = min(begindatetime))

dta_1 %>% 
ggplot(aes(begindatetime, yearday)) +
  geom_jitter() +
  annotate("text", label = "February", x = start$min, y = 50, size = 3, color = "black") +
  annotate("text", label = "March", x = start$min, y = 62, size = 3, color = "black") +
  annotate("text", label = "April", x = start$min, y = 104, size = 3, color = "black") +
  annotate("text", label = "May", x = start$min, y = 135, size = 3, color = "black") +
  annotate("text", label = "June", x = start$min, y = 155, size = 3, color = "black") +
  geom_hline(yintercept = c(60, 91, 121, 152), linetype = 2, color = "blue") +
  theme_bw() +
  scale_x_datetime() +
  ylim(45, 165) +
  labs(x = "Year",
       y = "Yearday",
       caption = "Blue dotted lines demarcate months.\n Leap year was not taken into account when determining yearday breaks for months.")


Questions:

  • Should the February and June data be removed?
  • Should I only be looking at the data from April?

Summary of SPECIES habitat

To recap from EDA 3:

The SPECIES table in FishBase has information on the depth range and the preferred environment of each species. The following definitions are from the website.

Habitat: Indicates the particular environment preferred by the species, with the following choices (adapted from Holthus and Maragos 1995):

pelagic: occurring mainly in the water column between 0 and 200 m, not feeding on benthic organisms;

benthopelagic: living and/or feeding on or near the bottom, as well as in midwater, between 0 and 200 m;

demersal: living and/or feeding on or near the bottom, between 0 and 200 m;

reef-associated: living and/or feeding on or near reefs, between 0 and 200 m;

bathypelagic: occurring mainly in open water below 200 m, not feeding on benthic organisms;

bathydemersal: living and/or feeding on or near the bottom, below 200 m.

Holthus, P.F. and J.E. Maragos. 1995. Marine ecosystem classification for the tropical island Pacific, p. 239-278. In J.E. Maragos, M.N.A. Peterson, L.G. Eldredge, J.E. Bardach and H.F. Takeuchi (eds.) Marine and coastal biodiversity in the tropical island Pacific region. Vol. 1. Species Management and Information Management Priorities. East-West Center, Honolulu, Hawaii. 424 p. 


Pelagic-neritic, pleagic-oceanic, and reef-associated are also choices for preferred depth but definitions are not given on the FishBase SPECIES table website.

Pelagic-neritic: the marine waters column above the continental shelf, some definitions specify the low tide line to the shelf break;

pleagic-oceanic: living in all open water that is not pelagic-neritic;

reef-associated: associated with reefs.

Counts for species

(Repeat from EDA 3)

These are the counts of species listed in the BTS data broken down by their habitat type and grouped by “fish” or “not fish”. These are not individual counts from the trawl data.

# setting up levels of habitat variables for plotting order
# hab_levs_1 <- c("reef-associated", "pelagic", "pelagic-neritic","pelagic-oceanic", 
#             "benthic", "benthopelagic", "demersal", "bathypelagic", "bathydemersal")

hab_levs_2 <- c("bathydemersal", "bathypelagic",  "demersal", "benthopelagic",
                "benthic", "pelagic-oceanic", "pelagic-neritic", "pelagic",
                "reef-associated")
labels <- c(fish = "Fish", not_fish = "Not Fish")

# plot of habitat types represented in the BTS data (species only)
sp_fish %>% 
  bind_rows(sp_sealife) %>% 
  count(type, DemersPelag) %>% 
  ggplot(aes(fct_relevel(DemersPelag, hab_levs_2), n, fill = type, label = n)) +
  geom_col() +
  geom_text(hjust = -0.5) +
  coord_flip() +
  facet_grid(~ type, labeller=labeller(type = labels)) +
  ylim(0, 150) +
  theme_bw() +
  theme(legend.position = "none") +
  scale_fill_viridis_d() +
  theme(axis.text.y = element_text(size = 10)) +
  labs(y = "Count",
       x = NULL, 
       caption = "Count of habitat type associated with species observations in the BTS (2008 - 2019) data")

The section below is code to show the counts from the graph above in tables. It is commented out. Delete the # in the code below if you would like to see the tables along with the graph.

# sp_fish %>% 
#   count(DemersPelag) %>% 
#   kable(col.names = c("Habitat", "Count"), 
#         caption = "Counts of species....") %>% 
#   kable_styling(bootstrap_options = c("striped", "hover"))
# 
# sp_sealife %>% 
#   count(DemersPelag) %>% 
#   kable(col.names = c("Habitat", "Count"), 
#         caption = "Counts of species....") %>% 
#   kable_styling(bootstrap_options = c("striped", "hover"))


Joining the trawl data to the species data

All the work so far has been to assess the species in the data sets without looking at the counts of organisms from the trawls. Here we add the taxonomic and habitat data to the trawl count data for species.

# making a single data frame to add to the count data
both_fish_sealife <- bind_rows(sp_fish, sp_sealife)

# combining the count data with the species data
dta_sp_only <- dta_ed %>% 
  left_join(both_fish_sealife, by = c("taxa"= "BTS_sp")) %>% # using the BTS names here
  rename(BTS_sp = taxa) %>% 
  filter(!is.na(Species)) %>% 
  mutate(level = "Species")

Note: dta_sp_only from EDA 2 was saved in the input_output folder as dta_and_species.rds. That file version does not have the swim bladder or habitat information from EDA 3 attached to it. Please be aware of what data you have loaded.

Broad scale counts

What are the broad-scale counts for “Fish” and “Not Fish” species by Cruise ID?

dta_sp_only %>% 
  count(type, cruiseid) %>% 
  ggplot(aes(type, n, fill = type, label = n)) + 
  geom_col(position =  position_dodge2(preserve = "single")) +
  geom_text(hjust = -0.1, size = 3) +
  coord_flip() +
  ylim(0, 7000) +
  scale_fill_viridis_d(breaks = c("not_fish", "fish"),
                       labels = c("Not Fish", "Fish")) +
  facet_wrap(~ cruiseid) +
  theme_bw() +
  theme(axis.text.y = element_blank(),
        axis.ticks.y = element_blank()) +
  labs(y = "Count",
       x = NULL,
       fill = NULL,
       caption = "Counts of fish individuals and not fish individuals from the BTS 2008 - 2019 data.")


Habitat types

What is the breakdown of habitat types for “Fish” and “Not Fish” species by Cruise ID?

dta_sp_only %>%
  count(DemersPelag, cruiseid, type) %>% 
  ggplot(aes(fct_relevel(DemersPelag, hab_levs_2), n, fill = type)) +
  geom_col(position =  position_dodge2(preserve = "single")) +
  scale_fill_viridis_d(labels = c("Fish", "Not Fish")) +
  coord_flip() +
  facet_wrap(~ cruiseid) +
  theme_bw() +
  labs(y = "Count",
       x = NULL,
       fill = NULL)

There are 338 species represented in the BTS data. Making a bar chart of all of them would be an exercise in very bad data visualization.

Question:

How would you like the data aggregated and visualized?

Guessing at groups


For the species that are Not Fish:

Vertebrates Invertebrates
Sea turtles Octopodes
Seals Squids
Crabs
Lobsters
Shrimps
Whelks
Scallop
Sea starts


For the species that are Fish:
- I need to look at the species list in EDA 2 and would like some help here please. :)

Quick and dirty example of guild plotting

I am pulling out everything in the trawl data from Class Elasmobranchii to play with.

sharks <- dta_sp_only %>% 
  filter(Class == "Elasmobranchii")

sharks %>% 
  count(BTS_sp, Com_name) %>% 
  kable(col.names = c("Scientific name", "Common Name", "Count")) %>% 
  kable_styling(bootstrap_options = c("striped", "hover")) %>% 
  column_spec(1, italic = TRUE)
Scientific name Common Name Count
Alopias superciliosus Bigeye thresher 1
Amblyraja radiata Starry ray 556
Carcharhinus acronotus Blacknose shark 1
Carcharhinus brevipinna Spinner shark 1
Carcharhinus falciformis Silky shark 2
Carcharhinus obscurus Dusky shark 4
Carcharhinus plumbeus Sandbar shark 16
Carcharias taurus Sand tiger shark 11
Cetorhinus maximus Basking shark 2
Dasyatis americana Southern stingray 14
Dasyatis centroura Roughtail stingray 57
Dasyatis say Bluntnose stingray 26
Dipturus laevis Barndoor skate 1003
Etmopterus gracilispinis Broadbanded lanternshark 3
Galeocerdo cuvier Tiger shark 1
Gymnura altavela Spiny butterfly ray 32
Gymnura micrura Smooth butterfly ray 14
Leucoraja erinacea Little skate 2391
Leucoraja garmani Rosette skate 370
Leucoraja ocellata Winter skate 2028
Malacoraja senta Smooth skate 748
Mustelus canis Dusky smooth-hound 356
Myliobatis freminvillei Bullnose eagle ray 15
Myliobatis goodei Southern eagle ray 2
Raja eglanteria Clearnose skate 543
Rhinoptera bonasus Cownose ray 1
Rhizoprionodon terraenovae Atlantic sharpnose shark 21
Scyliorhinus retifer Chain catshark 541
Sphyrna lewini Scalloped hammerhead 2
Sphyrna zygaena Smooth hammerhead 1
Squalus acanthias Picked dogfish 2419
Squatina dumeril Atlantic angel shark 131
Torpedo nobiliana Electric ray 60 
sharks %>% 
  count(Com_name, cruiseid) %>% 
  ggplot(aes(fct_reorder(Com_name, n), n)) +
  geom_col() +
  coord_flip() +
  facet_wrap(~ cruiseid) +
  theme_bw() +
  theme(axis.text.y = element_text(size = 7)) +
  labs(y = "Count", 
       x = NULL,
       caption = "Counts of Elasmobranchii catches by cruise.")

sharks %>%
  count(Com_name, month) %>%
  filter(!is.na(month)) %>% 
  ggplot(aes(fct_reorder(Com_name, n), n)) +
  geom_col() +
  coord_flip() +
  facet_wrap(~ month, scales = "free_x") +
  theme_bw() +
  theme(axis.text.y = element_text(size = 8)) +
  labs(y = "Count", 
       x = NULL,
       caption = "Counts of Elasmobranchii catches by month. \n x-axis on different scales.")

sharks %>% 
  count(Com_name, DemersPelag) %>% 
  filter(n > 250) %>% 
  ggplot(aes(fct_reorder(Com_name, n), n)) +
  geom_col() +
  facet_wrap(~ DemersPelag) +
  coord_flip() +
  theme_bw() +
  theme(axis.text.y = element_text(size = 8)) +
  labs(y = "Count", 
       x = NULL,
       fill = "Habitat",
       caption = "Counts greater than 250 individual Elasmobranchii.")

sharks %>% 
  count(Com_name, DemersPelag) %>% 
  filter(n <= 250) %>% 
  ggplot(aes(fct_reorder(Com_name, n), n)) +
  geom_col() +
  facet_wrap(~ DemersPelag) +
  coord_flip() +
  theme_bw() +
  theme(axis.text.y = element_text(size = 8)) +
  labs(y = "Count", 
       x = NULL,
       fill = "Habitat",
       caption = "Counts less than or equal to 250 individual Elasmobranchii.")

sharks %>% 
  filter(!is.na(year)) %>%
  mutate(year = as.factor(year)) %>% 
  count(Com_name, year) %>% 
  ggplot(aes(year, n, group = Com_name, color = Com_name)) +
  geom_line() +
  scale_color_viridis_d() +
  theme_bw() +
  theme(legend.position = "bottom") +
  guides(color = guide_legend(ncol = 4)) +
  labs(x = "Year",
       y = "Count",
       color = NULL,
       caption = "Amount of individuals collected during the trawl surveys for that year.")